void write_smOpt(const char *fname, int minia, int maxia){
	
	FILE * out_smmfile; out_smmfile  = fopen (fname, "w"); 
	    
    for (int ia = minia; ia <= maxia; ia++){ // NUM_age
        
        
        for (int iper = 0; iper < NUM_smmPerson; iper++){
                
		
		double theta[] = {sm_theta[0][iper], sm_theta[1][iper]}; 
			
		struct_infoset infoset;
		
		infoset.theta = theta; 
		infoset.health   = sm_cprice[ia][iper]; 
		infoset.logcprice   = log(sm_cprice[ia][iper]); 
				
		infoset.ia       = ia; 
		infoset.de_prev  = sm_de_prev[ia][iper];   
		infoset.debtlimitLag = 0.0; 
		if (ia>0) infoset.debtlimitLag = sm_debtlimit_p[ia-1][iper];
		
		infoset.educ     = sm_educ[ia][iper]; 
		infoset.kstock   = sm_kstock[ia][iper]; 
		infoset.qstock   = sm_qstock[ia][iper]; 
		infoset.savings  = sm_savings[ia][iper]; 
		
		infoset.ieduccatpr  = sm_ieduccatpr[iper]; 		
		infoset.isavingspr  = sm_isavingspr[iper];  
		
		infoset.eps_ue 	        = sm_rng_eps[ia][iper][1] ; 
		infoset.eps_uq          = sm_rng_eps[ia][iper][2] ; 
		
		infoset.eps_logwf       = get_rng_epsw(infoset.educ, sm_rng_eps[ia][iper][6]); 
				
		infoset.eps_health_next = sm_rng_eps[ia][iper][7];
		
		#ifndef INFOSET_hnext
						
		infoset.eps_health_next = 0.0; 
						
		#endif
						
		infoset.eps_ptr   		= sm_rng_eps[ia][iper][8]; 	 
		infoset.eps_ptrdummy   	= sm_rng_eps[ia][iper][9]; 
		
		infoset.skill           = get_skill(sm_educ[ia][iper], theta, ia, infoset.eps_logwf); 

						
			infoset.logcprice   = log(sm_cprice[ia][iper]) ; 			
			infoset.health   	= sm_cprice[ia][iper]; 
						
			infoset.savings = sm_savings[ia][iper]; 
        
		
			if (sm_isavingspr[iper] != 0) {printf("\n Error write_all.cpp 71 exiting"); exit(1); }		
 			const double delta   = get_delta(infoset.theta, infoset.qstock); 	
 			infoset.delta  = delta; 
 				
			double _util[2], _v[2], _evnext[2], _delta[2], _probs[2]; 
			for (int dq=0; dq<2; dq++){
				
				_v[dq] = get_valueEGM(sm_savings[ia+1][iper], sm_c[ia][iper], sm_de[ia][iper], sm_dk[ia][iper], dq, infoset, _util[dq], _evnext[dq], _probs[dq], _delta[dq]); 
					
			}
			
			double elife = get_lifeExpectancy_dynamic(ia, iper); 
						
			#ifdef USE_omp
            
			#pragma omp critical			
            
            #endif
            {
			
			fprintf(out_smmfile, "%f \t %f \t", elife, sm_probs[ia][iper]); 

		    fprintf(out_smmfile, "%d \t %d \t %f \t %f \t %.8f \t %2d \t %2f \t %.8f  \t", iper, ia, sm_theta[0][iper], sm_theta[1][iper], sm_cprice[ia][iper], sm_educ[ia][iper], sm_expr[ia][iper], sm_savings[ia][iper]);  // sm_c[ia][iper], 

		    fprintf(out_smmfile, "%d \t %f \t %d \t %.8f \t %.8f \t ", sm_de[ia][iper], sm_dk[ia][iper], sm_dq[ia][iper], sm_wage[ia][iper], sm_savings[ia+1][iper]) ; 
		    
		    fprintf(out_smmfile, "%d \t %d \t %f \t %f \t %f \t %f \t ", sm_ieduccatpr[iper], sm_add[ia][iper], sm_measure_thetac[iper], sm_measure_thetan[iper], sm_tax_cprice[ia][iper], sm_skill[ia][iper]) ; 
			
			fprintf(out_smmfile, "%f \t ", sm_earnings[ia][iper]) ; 
						
			fprintf(out_smmfile, "%f \t ", delta) ; 
 			
		    fprintf(out_smmfile, "%f \t ", sm_qstock[ia][iper]) ; 			
			
			fprintf(out_smmfile, "%f \t %f \t ", sm_value[ia][iper], infoset.eps_uq * sigeps_uq) ;
			for (int dq=0; dq<2; dq++){ fprintf(out_smmfile, "%f \t %f \t %f \t %f \t %f \t ", _v[dq], _util[dq], _evnext[dq], _delta[dq], _probs[dq]); }

            fprintf(out_smmfile, "%f \n", sm_wgt[iper]) ;
			
			}
				
		}
        
#ifdef USE_omp
        
#pragma omp barrier
#pragma omp flush
        
#endif
        
        
	}
	
	fclose(out_smmfile); 

}




int get_stat(const int i_n){
	
	int iage30 = 30 - INIT_age; 

	double elife = 0.0; 
	for (int iper = 0; iper < NUM_smmPerson; iper++){
		elife += get_lifeExpectancy_dynamic(0, iper); 
	}
	elife /= NUM_smmPerson; 
	
	AVEagedeath_a15 [i_n] = elife; 
	AVEyearsleft_a15[i_n] = elife - 1.0 - 15.0; 
	

	AVEdqever_a25[i_n] = my_countMore(NUM_smmPerson, sm_add[25-INIT_age], 0.5) / (1.0*NUM_smmPerson); 

	AVEdq_a15[i_n]		 = my_mean(NUM_smmPerson, sm_dq[0]);	
	AVEdq_a20[i_n]		 = my_mean(NUM_smmPerson, sm_dq[20-INIT_age]);
	AVEdq_a25[i_n]		 = my_mean(NUM_smmPerson, sm_dq[25-INIT_age]);
	AVEdq_a30[i_n]		 = my_mean(NUM_smmPerson, sm_dq[30-INIT_age]);	


	AVEhgc_a15[i_n]		 = my_mean(NUM_smmPerson, sm_educ[0]);	
	AVEhgc_a20[i_n]		 = my_mean(NUM_smmPerson, sm_educ[20-INIT_age]);
	AVEhgc_a25[i_n]		 = my_mean(NUM_smmPerson, sm_educ[25-INIT_age]);
	AVEhgc_a30[i_n]		 = my_mean(NUM_smmPerson, sm_educ[30-INIT_age]);	
	

	AVEc_a15[i_n]		 = my_mean(NUM_smmPerson, sm_c[0]);	
	AVEc_a20[i_n]		 = my_mean(NUM_smmPerson, sm_c[20-INIT_age]);
	AVEc_a25[i_n]		 = my_mean(NUM_smmPerson, sm_c[25-INIT_age]);
	AVEc_a30[i_n]		 = my_mean(NUM_smmPerson, sm_c[30-INIT_age]);	


	AVEnetworth_a15[i_n]		 = my_mean(NUM_smmPerson, sm_savings[0]);	
	AVEnetworth_a20[i_n]		 = my_mean(NUM_smmPerson, sm_savings[20-INIT_age]);
	AVEnetworth_a25[i_n]		 = my_mean(NUM_smmPerson, sm_savings[25-INIT_age]);
	AVEnetworth_a30[i_n]		 = my_mean(NUM_smmPerson, sm_savings[30-INIT_age]);	

	AVEearnings_a15[i_n]		 = my_mean(NUM_smmPerson, sm_earnings[0]);	
	AVEearnings_a20[i_n]		 = my_mean(NUM_smmPerson, sm_earnings[20-INIT_age]);
	AVEearnings_a25[i_n]		 = my_mean(NUM_smmPerson, sm_earnings[25-INIT_age]);
	AVEearnings_a30[i_n]		 = my_mean(NUM_smmPerson, sm_earnings[30-INIT_age]);	


	AVEvalue_a15[i_n]		 = my_mean(NUM_smmPerson, sm_value[0]);	
	AVEvalue_a20[i_n]		 = my_mean(NUM_smmPerson, sm_value[20-INIT_age]);
	AVEvalue_a25[i_n]		 = my_mean(NUM_smmPerson, sm_value[25-INIT_age]);
	AVEvalue_a30[i_n]		 = my_mean(NUM_smmPerson, sm_value[30-INIT_age]);	


	AVEqstock_a15[i_n]		 = my_mean(NUM_smmPerson, sm_qstock[0]);	
	AVEqstock_a20[i_n]		 = my_mean(NUM_smmPerson, sm_qstock[20-INIT_age]);
	AVEqstock_a25[i_n]		 = my_mean(NUM_smmPerson, sm_qstock[25-INIT_age]);
	AVEqstock_a30[i_n]		 = my_mean(NUM_smmPerson, sm_qstock[30-INIT_age]);
	
	
	AVEadd_a15[i_n]		 = my_mean(NUM_smmPerson, sm_add[0]);	
	AVEadd_a20[i_n]		 = my_mean(NUM_smmPerson, sm_add[20-INIT_age]);
	AVEadd_a25[i_n]		 = my_mean(NUM_smmPerson, sm_add[25-INIT_age]);
	AVEadd_a30[i_n]		 = my_mean(NUM_smmPerson, sm_add[30-INIT_age]);
		
		
	// ... start calculating statistics ... 	
		
	for (int i=0; i<4; i++) AVEeduccat4[i_n][i] = 0.0; 
	
	AVEeverdq[i_n] = 0.0; AVEinitdq[i_n] = 0.0; double n_initdq = 0.0; 	
	for (int iper = 0; iper < NUM_smmPerson; iper++){ 
		
		AVEeduccat4[i_n][0] += (sm_educ[30-INIT_age][iper]<= 11) / (1.0*NUM_smmPerson); 										 
		AVEeduccat4[i_n][1] += (sm_educ[30-INIT_age][iper] == 12) / (1.0*NUM_smmPerson); 										 
		AVEeduccat4[i_n][2] += (sm_educ[30-INIT_age][iper] >= 13) * (sm_educ[30-INIT_age][iper]<= 15) / (1.0*NUM_smmPerson); 	 
		AVEeduccat4[i_n][3] += (sm_educ[30-INIT_age][iper] >= 16) / (1.0*NUM_smmPerson); 										 
				
		
		AVEeverdq[i_n] += ( sm_add[iage30][iper] > 0.5 ) / (1.0*NUM_smmPerson);

		AVEinitdq[i_n] += 1.0 * sm_dq[0][iper] * (sm_add[0][iper] < 0.5); 
		n_initdq += 1.0 * (sm_add[0][iper] < 0.5); 		
	}
	AVEinitdq[i_n] = AVEinitdq[i_n]/n_initdq; 
//	printf("\n check0: %f ", AVEinitdq[i_n]);
	
	AVEinitdq[i_n] = my_mean_conditionalLess(NUM_smmPerson, sm_dq[0], sm_add[0], 0.01);
//	printf("\n check1: %f ", AVEinitdq[i_n]); 	

	
	AVEdq1630[i_n] = 0.0; 
	for (int ia=0; ia<=iage30; ia++){
		AVEdq1630[i_n] += my_mean(NUM_smmPerson, sm_dq[ia]) / (1.0*(iage30+1.0));	
	}
	
			
	Probhgc12_a19[i_n] = my_countMore(NUM_smmPerson, sm_educ[19-INIT_age], 11) / (1.0*NUM_smmPerson) ; 
	Probhgc13_a21[i_n] = my_countMore(NUM_smmPerson, sm_educ[21-INIT_age], 12) / (1.0*NUM_smmPerson) ; 	
	Probhgc16_a25[i_n] = my_countMore(NUM_smmPerson, sm_educ[25-INIT_age], 15) / (1.0*NUM_smmPerson) ; 
	
	AVEhgc[i_n]		 = my_mean(NUM_smmPerson, sm_educ[iage30]);
	AVEhealth[i_n] 	 = my_mean(NUM_smmPerson, sm_cprice[iage30]);
	
	AVEskill[i_n]    = my_mean(NUM_smmPerson, sm_skill[iage30]);
	AVEqstock[i_n]   = my_mean(NUM_smmPerson, sm_qstock[iage30]);
	AVEdq[i_n]		 = my_mean(NUM_smmPerson, sm_dq[iage30]);	
	
	
	AVEearnings[i_n] = my_mean(NUM_smmPerson, sm_earnings[iage30]);
	AVEc[i_n] 		 = my_mean(NUM_smmPerson, sm_c[iage30]);
	AVEnetworth[i_n] = my_mean(NUM_smmPerson, sm_savings[iage30]);
		
	
	AVEinitvalue[i_n]    = my_mean(NUM_smmPerson, sm_value[0]);    
	AVEvaluea30[i_n]    = my_mean(NUM_smmPerson, sm_value[iage30]); 	
	
	
	
	WorkpartEnroll[i_n]  = SMM_moments[13];	
	
	
	GiniWealth[i_n]      = my_gini(NUM_smmPerson, sm_savings[iage30]);     
    GiniEarnings[i_n]    = my_gini(NUM_smmPerson, sm_earnings[iage30]); 
	
	VARlogearnings[i_n] = my_varlog(NUM_smmPerson, sm_earnings[iage30]);  
	VARlogc[i_n]        = my_varlog(NUM_smmPerson, sm_c[iage30]); 	
	VARlogskill[i_n]    = my_varlog(NUM_smmPerson, sm_skill[iage30]); 		
	VARloghgc[i_n]      = my_varlog(NUM_smmPerson, sm_educ[iage30]); 
	VARlogcprice[i_n]    = my_varlog(NUM_smmPerson, sm_cprice[iage30]); 		
	VARlogqstock[i_n]    = my_varlog1(NUM_smmPerson, sm_qstock[iage30]); 
	VARloginitvalue[i_n]= my_varlog(NUM_smmPerson, sm_value[0]);  
	VARlogvaluea30[i_n]= my_varlog(NUM_smmPerson, sm_value[iage30]); 
	

	
	
	{
	int icat; int hgcmin, hgcmax; 
	
	icat = 1; 
	hgcmin = 0; hgcmax = 11; 
	
	educcatAVEinitvalue[i_n][icat-1]	   = my_mean_conditional(NUM_smmPerson, sm_value[0], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEvaluea30[i_n][icat-1]	    = my_mean_conditional(NUM_smmPerson, sm_value[iage30], sm_educ[iage30], hgcmin, hgcmax); 
	educcatAVEhealth[i_n][icat-1]	 	 = my_mean_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEqstock[i_n][icat-1]	   	 = my_mean_conditional(NUM_smmPerson, sm_qstock[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEdq[i_n][icat-1]			 = my_mean_conditional(NUM_smmPerson, sm_dq[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEearnings[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEc[i_n][icat-1]	 		 = my_mean_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEnetworth[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_savings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	
	educcatAVElogcprice[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogearnings[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogc[i_n][icat-1]	 		 = my_meanlog_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);

	
	icat = 2; 
	hgcmin = 12; hgcmax = 12; 
	
	educcatAVEinitvalue[i_n][icat-1]	   = my_mean_conditional(NUM_smmPerson, sm_value[0], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEvaluea30[i_n][icat-1]	    = my_mean_conditional(NUM_smmPerson, sm_value[iage30], sm_educ[iage30], hgcmin, hgcmax); 
	educcatAVEhealth[i_n][icat-1]	 	 = my_mean_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEqstock[i_n][icat-1]	   	 = my_mean_conditional(NUM_smmPerson, sm_qstock[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEdq[i_n][icat-1]			 = my_mean_conditional(NUM_smmPerson, sm_dq[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEearnings[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEc[i_n][icat-1]	 		 = my_mean_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEnetworth[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_savings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	
	educcatAVElogcprice[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogearnings[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogc[i_n][icat-1]	 		 = my_meanlog_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);
	
	
	icat = 3; 
	hgcmin = 13; hgcmax = 15; 
	
	educcatAVEinitvalue[i_n][icat-1]	   = my_mean_conditional(NUM_smmPerson, sm_value[0], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEvaluea30[i_n][icat-1]	    = my_mean_conditional(NUM_smmPerson, sm_value[iage30], sm_educ[iage30], hgcmin, hgcmax); 
	educcatAVEhealth[i_n][icat-1]	 	 = my_mean_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEqstock[i_n][icat-1]	   	 = my_mean_conditional(NUM_smmPerson, sm_qstock[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEdq[i_n][icat-1]			 = my_mean_conditional(NUM_smmPerson, sm_dq[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEearnings[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEc[i_n][icat-1]	 		 = my_mean_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEnetworth[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_savings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	
	educcatAVElogcprice[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogearnings[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogc[i_n][icat-1]	 		 = my_meanlog_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);


	icat = 4; 
	hgcmin = 16; hgcmax = 20; 
	
	educcatAVEinitvalue[i_n][icat-1]	   = my_mean_conditional(NUM_smmPerson, sm_value[0], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEvaluea30[i_n][icat-1]	    = my_mean_conditional(NUM_smmPerson, sm_value[iage30], sm_educ[iage30], hgcmin, hgcmax); 
	educcatAVEhealth[i_n][icat-1]	 	 = my_mean_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEqstock[i_n][icat-1]	   	 = my_mean_conditional(NUM_smmPerson, sm_qstock[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEdq[i_n][icat-1]			 = my_mean_conditional(NUM_smmPerson, sm_dq[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEearnings[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEc[i_n][icat-1]	 		 = my_mean_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVEnetworth[i_n][icat-1]	 = my_mean_conditional(NUM_smmPerson, sm_savings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	
	educcatAVElogcprice[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_cprice[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogearnings[i_n][icat-1]	 = my_meanlog_conditional(NUM_smmPerson, sm_earnings[iage30], sm_educ[iage30], hgcmin, hgcmax);
	educcatAVElogc[i_n][icat-1]	 		 = my_meanlog_conditional(NUM_smmPerson, sm_c[iage30], sm_educ[iage30], hgcmin, hgcmax);

 	}

	
	// ... STATISTICS AMONG THOSE WHO WERE EVER SMOKERS in the benchmark model 	
	
	int dummy_everdq[NUM_smmPerson]; int dummy_dqa30dq[NUM_smmPerson]; 
		
	if (i_n == 0){ // write benchmark ever smokers 
	
		// only update at the benchmark model 
		for (int iper = 0; iper < NUM_smmPerson; iper++){
		dummy_everdq[iper] = ( sm_add[30-INIT_age][iper] + sm_dq[30-INIT_age][iper] > 0.5 ); 
		}
		printvec("out_files/dummy_everdq.txt", dummy_everdq, NUM_smmPerson); 
		printvec("out_files/addiction_a30.txt", sm_add[30-INIT_age], NUM_smmPerson); 
		printvec("out_files/addiction_a25.txt", sm_add[25-INIT_age], NUM_smmPerson); 		
		printvec("out_files/dummy_dqa30dq.txt", sm_dq[30-INIT_age], NUM_smmPerson); 
		printvec("out_files/dummy_dqa25dq.txt", sm_dq[25-INIT_age], NUM_smmPerson); 
				
	} else { // read benchmark ever smokers 
	
		// ... read dummy_everdq from benchmark model 
		ifstream infile; 
	    infile.open("out_files/dummy_everdq.txt"); 	
	    if(!infile) {cerr<<"\n Error reading the file out_files/dummy_everdq.txt \n"; exit(1); }
		for (int iper = 0; iper < NUM_smmPerson; iper++){
		infile >> dummy_everdq[iper]; 
		}
		infile.close(); 

	    infile.open("out_files/dummy_dqa25dq.txt"); 	
	    if(!infile) {cerr<<"\n Error reading the file out_files/dummy_dqa25dq.txt \n"; exit(1); }
		for (int iper = 0; iper < NUM_smmPerson; iper++){
		infile >> dummy_dqa30dq[iper]; 
		}
		infile.close(); 	
	}
	
	
	double ever_elife = 0.0; double n_ever = 0.0; 
	for (int iper = 0; iper < NUM_smmPerson; iper++){
		if (dummy_everdq[iper]==1){

		ever_elife += get_lifeExpectancy_dynamic(0, iper); 
		n_ever     += 1.0; 
			
		}
	}
	ever_elife /= n_ever ; 
	
	everAVEagedeath_a15 [i_n] = ever_elife; 
	everAVEyearsleft_a15[i_n] = ever_elife - 1.0 - 15.0; 
	
	
	
	
  int i2530 = 1; 
	
	AVEhgc11less [i_n]= my_mean(NUM_smmPerson, sm_educcat[i2530][0]);            
	AVEhgc12     [i_n]= my_mean(NUM_smmPerson, sm_educcat[i2530][1]); 
	AVEhgc1315   [i_n]= my_mean(NUM_smmPerson, sm_educcat[i2530][2]);             
	AVEhgc16more [i_n]= my_mean(NUM_smmPerson, sm_educcat[i2530][3]);

	AVEhgc11less_everdq [i_n]= my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][0], dummy_everdq, 1,  1);        
	AVEhgc12_everdq     [i_n]= my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][1], dummy_everdq, 1,  1); 
	AVEhgc1315_everdq   [i_n]= my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][2], dummy_everdq, 1,  1);           
	AVEhgc16more_everdq [i_n]= my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][3], dummy_everdq, 1,  1); 
	
	AVEhgc12more        [i_n]= 1.0 - AVEhgc11less[i_n];            
	AVEhgc12more_everdq [i_n]= 1.0 - AVEhgc11less_everdq[i_n]; 
	
	iage30 = 30 - INIT_age; 
			
	const int sum_dummy_everdq = my_sum(NUM_smmPerson, dummy_everdq); 

	everAVEdq_a15[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_dq[15-INIT_age], dummy_everdq, 0.9);
	everAVEdq_a20[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_dq[20-INIT_age], dummy_everdq, 0.9);
	everAVEdq_a25[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_dq[25-INIT_age], dummy_everdq, 0.9);
	everAVEdq_a30[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_dq[30-INIT_age], dummy_everdq, 0.9);


	everAVEadd[i_n] 	 = my_mean_conditionalMore(NUM_smmPerson, sm_add[iage30], dummy_everdq, 0.9);
	everAVEqstock[i_n]   = my_mean_conditionalMore(NUM_smmPerson, sm_qstock[iage30], dummy_everdq, 0.9);
	everAVEdq[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_dq[iage30], dummy_everdq, 0.9);	

	everAVEhgc[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_educ[iage30], dummy_everdq, 0.9);
	everAVEearnings[i_n] = my_mean_conditionalMore(NUM_smmPerson, sm_earnings[iage30], dummy_everdq, 0.9);
	everAVEc[i_n] 		 = my_mean_conditionalMore(NUM_smmPerson, sm_c[iage30], dummy_everdq, 0.9);
	everAVEnetworth[i_n] = my_mean_conditionalMore(NUM_smmPerson, sm_savings[iage30], dummy_everdq, 0.9);
		
		
	everAVEeverdq[i_n] = 0.0; 
	everProbhgc12_a19[i_n] = 0.0; 
	everProbhgc13_a21[i_n] = 0.0; 
	everProbhgc16_a25[i_n] = 0.0; 
	for (int iper = 0; iper < NUM_smmPerson; iper++){
		
		
		everProbhgc12_a19[i_n] += (sm_educ[19-INIT_age][iper] > 11) * (dummy_everdq[iper]==1) / (1.0*sum_dummy_everdq); 
		everProbhgc13_a21[i_n] += (sm_educ[21-INIT_age][iper] > 12) * (dummy_everdq[iper]==1) / (1.0*sum_dummy_everdq); 
		everProbhgc16_a25[i_n] += (sm_educ[25-INIT_age][iper] > 15) * (dummy_everdq[iper]==1) / (1.0*sum_dummy_everdq); 

		everAVEeverdq[i_n] += ( sm_add[iage30][iper] > 0.5 ) * (dummy_everdq[iper]==1) / (1.0*sum_dummy_everdq); 
	
	}
	
	everAVEdq1630[i_n] = 0.0; 
	for (int ia=0; ia<=iage30; ia++){
		everAVEdq1630[i_n] += my_mean_conditionalMore(NUM_smmPerson, sm_dq[ia], dummy_everdq, 0.9) / (1.0*(iage30+1.0));	
	}
	
	everAVEskill[i_n]    = my_mean_conditionalMore(NUM_smmPerson, sm_skill[iage30], dummy_everdq, 0.9);
	everAVEinitvalue[i_n]   = my_mean_conditionalMore(NUM_smmPerson, sm_value[0], dummy_everdq, 0.9);    
	everAVEvaluea30[i_n]    = my_mean_conditionalMore(NUM_smmPerson, sm_value[iage30], dummy_everdq, 0.9);  	



	const int sum_dummy_dqa30dq = my_sum(NUM_smmPerson, dummy_dqa30dq); 

	dqa30AVEadd[i_n] 	 = my_mean_conditionalMore(NUM_smmPerson, sm_add[iage30], dummy_dqa30dq, 0.9);
	dqa30AVEqstock[i_n]   = my_mean_conditionalMore(NUM_smmPerson, sm_qstock[iage30], dummy_dqa30dq, 0.9);
	dqa30AVEdq[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_dq[iage30], dummy_dqa30dq, 0.9);	

	dqa30AVEhgc[i_n]		 = my_mean_conditionalMore(NUM_smmPerson, sm_educ[iage30], dummy_dqa30dq, 0.9);
	dqa30AVEearnings[i_n] = my_mean_conditionalMore(NUM_smmPerson, sm_earnings[iage30], dummy_dqa30dq, 0.9);
	dqa30AVEc[i_n] 		 = my_mean_conditionalMore(NUM_smmPerson, sm_c[iage30], dummy_dqa30dq, 0.9);
	dqa30AVEnetworth[i_n] = my_mean_conditionalMore(NUM_smmPerson, sm_savings[iage30], dummy_dqa30dq, 0.9);
		
		
	dqa30AVEdqa30dq[i_n] = 0.0; 
	dqa30Probhgc12_a19[i_n] = 0.0; 
	dqa30Probhgc13_a21[i_n] = 0.0; 
	dqa30Probhgc16_a25[i_n] = 0.0; 
	for (int iper = 0; iper < NUM_smmPerson; iper++){
		
		
		dqa30Probhgc12_a19[i_n] += (sm_educ[19-INIT_age][iper] > 11) * (dummy_dqa30dq[iper]==1) / (1.0*sum_dummy_dqa30dq); 
		dqa30Probhgc13_a21[i_n] += (sm_educ[21-INIT_age][iper] > 12) * (dummy_dqa30dq[iper]==1) / (1.0*sum_dummy_dqa30dq); 
		dqa30Probhgc16_a25[i_n] += (sm_educ[25-INIT_age][iper] > 15) * (dummy_dqa30dq[iper]==1) / (1.0*sum_dummy_dqa30dq); 

		dqa30AVEdqa30dq[i_n] += ( sm_dq[30-INIT_age][iper] ) * (dummy_dqa30dq[iper]==1) / (1.0*sum_dummy_dqa30dq); 
	
	}
	
	dqa30AVEdq1630[i_n] = 0.0; 
	for (int ia=0; ia<=iage30; ia++){
		dqa30AVEdq1630[i_n] += my_mean_conditionalMore(NUM_smmPerson, sm_dq[ia], dummy_dqa30dq, 0.9) / (1.0*(iage30+1.0));	
	}
	
	dqa30AVEskill[i_n]    = my_mean_conditionalMore(NUM_smmPerson, sm_skill[iage30], dummy_dqa30dq, 0.9);
	dqa30AVEinitvalue[i_n]   = my_mean_conditionalMore(NUM_smmPerson, sm_value[0], dummy_dqa30dq, 0.9);    
	dqa30AVEvaluea30[i_n]    = my_mean_conditionalMore(NUM_smmPerson, sm_value[iage30], dummy_dqa30dq, 0.9);  	
	
	iage30 = 30 - INIT_age; 
	
	
	
	return 0; 
}



int write_input(int parafit){ 
// ... parafit == 1: write parameters and fitted moments; 



	string vec_strNumb[] = { "CFZero", "CFOne", "CFTwo", "CFThree", "CFFour", "CFFive", "CFSix", "CFSeven", "CFEight", "CFNine"}; 		
   	string str; 

	int ncf = sizeof(vec_strNumb)/sizeof(string); 
    if (Num_cf != ncf) {printf("\n\n Errors in write_input: Num_cf %d, ncf %d ... exiting ...\n\n", Num_cf, ncf); exit(1); }
		
		
	
	int i = 0; 
	
	if (parafit == 1){	
	// ... Writing the output ... 
		
	
	printf("\n *** Writing parameters and model fit ***\n");
	
	FILE * input_estpar; 
	input_estpar  = fopen ("input_estpar.tex", "w"); 
	
	double sdthetac = sqrt( my_var(NUM_smmPerson, sm_theta[0]) ); 
	double sdthetan = sqrt( my_var(NUM_smmPerson, sm_theta[1]) );

	fprintf(input_estpar, "	\\newcommand{		\\sdthetac	}{  %.2f } \n",		sdthetac	 );
	fprintf(input_estpar, "	\\newcommand{		\\sdthetan	}{  %.2f } \n",		sdthetan	 );
	
	double newc[NUM_smmPerson]; int idq=0; 
	for (int iper = 0; iper < NUM_smmPerson; iper++){ 
		if (sm_dq[18-INIT_age][iper]==1) {
			newc[idq] = sm_c[18-INIT_age][iper]; 
			idq++; 
		}
	}
	double medconsp_smoker_a18 = my_median(idq, newc);
	double meanconsp_smoker_a18 = my_mean(idq, newc);	
	fprintf(input_estpar, "	\\newcommand{		\\MedConsumptionSmoker	}{  %.3f } \n",	medconsp_smoker_a18/1000.0 );
	fprintf(input_estpar, "	\\newcommand{		\\MeanConsumptionSmoker	}{  %.3f } \n",	meanconsp_smoker_a18/1000.0 );	
	
	
	fprintf(input_estpar, "	\\newcommand{		\\PctEffectclogwalphacCLG	}{  %.1f } \n",		100*(exp(clogw_alpha[0][3] * sdthetac) - 1.0) ); 
	fprintf(input_estpar, "	\\newcommand{		\\PctEffectclogwalphanCLG	}{  %.1f } \n",		100*(exp(clogw_alpha[1][3] * sdthetan) - 1.0) ); 
	fprintf(input_estpar, "	\\newcommand{		\\PctEffectclogwPbetaOne  	}{  %.1f } \n",		100*(exp(clogw_beta[1])-1.0)  	 ); 
	
	
	fprintf(input_estpar, "	\\newcommand{		\\PctRatiocuqalphaOneToOne	}{  %.0f } \n",		100*fabs(cu_q_alpha[1]*sdthetan/cu_q[1]) ); 
	fprintf(input_estpar, "	\\newcommand{		\\RatiocuqOneToe	}{  %.0f } \n",		        fabs(cu_q[1]/cu_qe_cat) ); 
	fprintf(input_estpar, "	\\newcommand{		\\RatiocuqOneTosq	}{  %.0f } \n",		        fabs(cu_q[1]/cu_qstock) ); 

	fprintf(input_estpar, "	\\newcommand{		\\RatiocuealphaZeroToOne	}{  %.0f } \n",		cu_e_alpha[0]*sdthetac/(cu_e_alpha[1]*sdthetan)	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\RatiocuealphaZeroToqstock	}{  %.0f } \n",		fabs(cu_e_alpha[0]*sdthetac/cu_e_qstock)	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\RatiocuePrToalphaZero	}{  %.0f } \n",		    cu_e_pr/(cu_e_alpha[0]*sdthetac)	 ); 
	
	
	fprintf(input_estpar, "	\\newcommand{		\\cuqOneStocksq	}{  %.4f } \n",		cu_q[1] + cu_qstock ); 	
	
	fprintf(input_estpar, "	\\newcommand{		\\cuqqpp  	}{  %.0f } \n",	cu_qq*100.0	); 	 
	
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaOneOneSD	}{  %.4f } \n",		delta1[1]*sdthetac	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaOneTwoSD	}{  %.4f } \n",		delta1[2]*sdthetan	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\PctcdeltaTwoZero	}{  %.3f } \n",		100.0*delta2[0]	 );	
	
	fprintf(input_estpar, "	\\newcommand{		\\cuealphaZeroSD	}{  %.4f } \n",		cu_e_alpha[0]*sdthetac	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuealphaOneSD 	}{  %.4f } \n",		cu_e_alpha[1]*sdthetan	 ); 	
	fprintf(input_estpar, "	\\newcommand{		\\cuqalphaZeroSD	}{  %.5f } \n",		cu_q_alpha[0]*sdthetac	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuqalphaOneSD 	}{  %.5f } \n",		cu_q_alpha[1]*sdthetan	 ); 	

	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacHSDSD 	}{  %.4f } \n",		clogw_alpha[0][0]*sdthetac	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanHSDSD 	}{  %.4f } \n",		clogw_alpha[1][0]*sdthetan	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacHSGSD	}{  %.4f } \n",		clogw_alpha[0][1]*sdthetac	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanHSGSD	}{  %.4f } \n",		clogw_alpha[1][1]*sdthetan	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacSCLSD	}{  %.4f } \n",		clogw_alpha[0][2]*sdthetac	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanSCLSD	}{  %.4f } \n",		clogw_alpha[1][2]*sdthetan	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacCLGSD	}{  %.4f } \n",		clogw_alpha[0][3]*sdthetac 	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanCLGSD	}{  %.4f } \n",		clogw_alpha[1][3]*sdthetan	 ); 


	fprintf(input_estpar, "	\\newcommand{		\\cudedq 	}{  %.4f } \n",		cu_de_dq 	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cueqstock 	}{  %.4f } \n",		cu_e_qstock 	 ); 
	
	fprintf(input_estpar, "	\\newcommand{		\\cug	}{  %.4f } \n",		cu_g	 );
	fprintf(input_estpar, "	\\newcommand{		\\clogwGraduate	}{  %.4f } \n",		clogw_skill_graduate	 );
	fprintf(input_estpar, "	\\newcommand{		\\clogwksqureGraduate	}{  %.4f } \n",		clogw_ksquare_graduate	 );

   
    
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaTwoZero	}{  %.5f } \n",		delta2[0]	 );
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaTwoOne	}{  %.4f } \n",		delta2[1]	 );
	
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaZero 	}{  %.4f } \n",		delta0	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaOneZero 	}{  %.4f } \n",		delta1[0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaOneOne	}{  %.4f } \n",		delta1[1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cdeltaOneTwo	}{  %.4f } \n",		delta1[2]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cumin 	}{  %.0f } \n",		cu_min 	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cueZero	}{  %.4f } \n",		cu_e[0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cueOne	}{  %.4f } \n",		cu_e[1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cueCLG	}{  %.4f } \n",		cu_e[2]	 );
    fprintf(input_estpar, " \\newcommand{       \\cueThree    }{  %.4f } \n",        cu_e[3]     );
	fprintf(input_estpar, "	\\newcommand{		\\cueReturn	}{  %.4f } \n",		cu_e_rcost	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuealphaZero	}{  %.4f } \n",		cu_e_alpha[0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuealphaOne 	}{  %.4f } \n",		cu_e_alpha[1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cueprZero	}{  %.4f } \n",		cu_e_pr	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cueage	}{  %.4f } \n",		cu_e_a	 ); 
	
        
	fprintf(input_estpar, "	\\newcommand{		\\sigepsue 	}{  %.4f } \n",		sigeps_ue	 ); 

    fprintf(input_estpar, " \\newcommand{        \\cukage    }{  %.4f } \n",        cu_ka      );

	fprintf(input_estpar, "	\\newcommand{		\\cukeHS	}{  %.4f } \n",		cu_ke[0]	 ); 	
	fprintf(input_estpar, "	\\newcommand{		\\cukeCG	}{  %.4f } \n",		cu_ke[1]	 ); 
	

	fprintf(input_estpar, "	\\newcommand{		\\cuqq	}{  %.4f } \n",		cu_qq	 );
	fprintf(input_estpar, "	\\newcommand{		\\cuqOne	}{  %.4f } \n",		cu_q[1]	 ); 
	
	
			
	fprintf(input_estpar, "	\\newcommand{		\\cuqZero	}{  %.4f } \n",		cu_q[0]	 ); 

	fprintf(input_estpar, "	\\newcommand{		\\cuqe	}{  %.4f } \n",		cu_qe_cat	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuqalphaZero	}{  %.5f } \n",		cu_q_alpha[0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuqalphaOne 	}{  %.5f } \n",		cu_q_alpha[1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\sigepsuq 	}{  %.4f } \n",		sigeps_uq	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\cuqstocksq  	}{  %.4f } \n",		cu_qstock 	 ); 
	
	
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaOne	}{  %.4f } \n",		clogw_skill[0]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaTwo  	}{  %.4f } \n",		clogw_skill[1]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaThree  	}{  %.4f } \n",		clogw_skill[2]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaFour  	}{  %.4f } \n",		clogw_skill[3]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaFive	}{  %.4f } \n",		clogw_skill[4]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaSix	}{  %.4f } \n",		clogw_skill[5]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwbetaSeven	}{  %.4f } \n",		clogw_skill[6]  	 ); 

	fprintf(input_estpar, "	\\newcommand{		\\clogwFbetaZero  	}{  %.4f } \n",		clogw_beta[0]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwPbetaOne  	}{  %.4f } \n",		clogw_beta[1]  	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwPbetaEnroll	}{  %.4f } \n",		clogw_beta[2] 	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwPbetaAge	}{  %.4f } \n",		clogw_beta[3]  	 ); 

	fprintf(input_estpar, "	\\newcommand{		\\sigepslogwNdk	}{  %.4f } \n",		clogw_ndk	 ); 

	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacHSD 	}{  %.4f } \n",		clogw_alpha[0][0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanHSD 	}{  %.4f } \n",		clogw_alpha[1][0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacHSG	}{  %.4f } \n",		clogw_alpha[0][1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanHSG	}{  %.4f } \n",		clogw_alpha[1][1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacSCL	}{  %.4f } \n",		clogw_alpha[0][2]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanSCL	}{  %.4f } \n",		clogw_alpha[1][2]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphacCLG	}{  %.4f } \n",		clogw_alpha[0][3] 	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\clogwalphanCLG	}{  %.4f } \n",		clogw_alpha[1][3]	 ); 



	
	
	fprintf(input_estpar, "	\\newcommand{		\\logwscaleHSD   	}{  %.4f } \n",		sigeps_logwscale[0]		 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwscaleHSG	}{  %.4f } \n",		sigeps_logwscale[1]    	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwscaleSCL	}{  %.4f } \n",		sigeps_logwscale[2]    	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwscaleCLG	}{  %.4f } \n",		sigeps_logwscale[3]    	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwscalePro	}{  %.4f } \n",		sigeps_logwscale[4]    	 ); 

	fprintf(input_estpar, "	\\newcommand{		\\logwshapeHSD	}{  %.4f } \n",		sigeps_logwshape[0]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwshapeHSG	}{  %.4f } \n",		sigeps_logwshape[1]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwshapeSCL	}{  %.4f } \n",		sigeps_logwshape[2]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwshapeCLG	}{  %.4f } \n",		sigeps_logwshape[3]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\logwshapePro	}{  %.4f } \n",		sigeps_logwshape[4]	 ); 
	
			
	fprintf(input_estpar, "	\\newcommand{		\\minepslogwHSD	}{  %.4f } \n",		MIN_epsw[0]   	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\minepslogwHSG	}{  %.4f } \n",		MIN_epsw[1]   	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\minepslogwSCL	}{  %.4f } \n",		MIN_epsw[2]   	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\minepslogwCLG	}{  %.4f } \n",		MIN_epsw[3]	 ); 
	fprintf(input_estpar, "	\\newcommand{		\\minepslogwPro	}{  %.4f } \n",		MIN_epsw[4]	 ); 
	

    double theta0[]={0,0};
	double cdeltaMean = get_delta(theta0, 0.0);
	fprintf(input_estpar, "	\\newcommand{	\\cexpdeltazero	}{  %.4f } \n", cdeltaMean );

	
	
	fclose(input_estpar); 
	

	
	FILE * input_vec_mean; 
	input_vec_mean  = fopen ("input_vec_mean.tex", "w"); 
	

	string count_str[] = { 
   "SMMzero	", 
   "SMMone	", 
   "SMMtwo	", 
   "SMMthree	", 
   "SMMfour	", 
   "SMMfive	", 
   "SMMsix	", 
   "SMMseven	", 
   "SMMeight	", 
   "SMMnine	", 
 
   "SMMten	", 
   "SMMeleven	", 
   "SMMtwelve	", 
   "SMMthirteen	", 
   "SMMfourteen	", 
   "SMMfifteen	", 
   "SMMsixteen	", 
   "SMMseventeen	", 
   "SMMeighteen	", 
   "SMMnineteen	", 
 
   "SMMtwenty	", 
   "SMMtwentyone	", 
   "SMMtwentytwo	", 
   "SMMtwentythree	", 
   "SMMtwentyfour	", 
   "SMMtwentyfive	", 
   "SMMtwentysix	", 
   "SMMtwentyseven	", 
   "SMMtwentyeight	", 
   "SMMtwentynine	", 
 
   "SMMthirty	", 
   "SMMthirtyone	", 
   "SMMthirtytwo	", 
   "SMMthirtythree	", 
   "SMMthirtyfour	", 
   "SMMthirtyfive	", 
   "SMMthirtysix	", 
   "SMMthirtyseven	", 
   "SMMthirtyeight	", 
   "SMMthirtynine	", 
 
   "SMMforty	", 
   "SMMfortyone	", 
   "SMMfortytwo	", 
   "SMMfortythree	", 
   "SMMfortyfour	", 
   "SMMfortyfive	", 
   "SMMfortysix	", 
   "SMMfortyseven	", 
   "SMMfortyeight	", 
   "SMMfortynine	", 
 
   "SMMfifty	", 
   "SMMfiftyone	", 
   "SMMfiftytwo	", 
   "SMMfiftythree	", 
   "SMMfiftyfour	", 
   "SMMfiftyfive	", 
   "SMMfiftysix	", 
   "SMMfiftyseven	", 
   "SMMfiftyeight	", 
   "SMMfiftynine	", 
 
   "SMMsixty	", 
   "SMMsixtyone	", 
   "SMMsixtytwo	", 
   "SMMsixtythree	", 
   "SMMsixtyfour	", 
   "SMMsixtyfive	", 
   "SMMsixtysix	", 
   "SMMsixtyseven	", 
   "SMMsixtyeight	", 
   "SMMsixtynine	", 
 
   "SMMseventy	", 
   "SMMseventyone	", 
   "SMMseventytwo	", 
   "SMMseventythree	", 
   "SMMseventyfour	", 
   "SMMseventyfive	", 
   "SMMseventysix	", 
   "SMMseventyseven	", 
   "SMMseventyeight	", 
   "SMMseventynine	", 
 
   "SMMeighty	", 
   "SMMeightyone	", 
   "SMMeightytwo	", 
   "SMMeightythree	", 
   "SMMeightyfour	", 
   "SMMeightyfive	", 
   "SMMeightysix	", 
   "SMMeightyseven	", 
   "SMMeightyeight	", 
   "SMMeightynine	", 
 
   "SMMninety	", 
   "SMMninetyone	", 
   "SMMninetytwo	", 
   "SMMninetythree	", 
   "SMMninetyfour	", 
   "SMMninetyfive	", 
   "SMMninetysix	", 
   "SMMninetyseven	", 
   "SMMninetyeight	", 
   "SMMninetynine	", 
   "SMMhundred	"  
	}; 

	int ncount = sizeof(count_str)/sizeof(string); 
	
	
	fprintf(input_vec_mean, "	\\newcommand{	\\SMMsumf 	}{  %.4f } \n",		SMM_sumf2	 ); 	i++; 

	i = 0; 
	
	for(int istat=0; istat < ncount; istat++){
	
	string str = count_str[istat]; 
	
	fprintf(input_vec_mean,  "	\\newcommand{	\\%s	}{  %.4f } \n",	str.c_str(), SMM_moments[i] ); i++; 
			
	}

	fclose(input_vec_mean); 	
	

	}  // parafit == 1
	else { 


	string vec_str[] = {
		
		"AVEagedeath", 
		"AVEyearsleft", 

		"everAVEagedeath", 
		"everAVEyearsleft", 
		
		"AVEhgcSCLmore", 
		"everAVEhgcSCLmore", 
				
		"AVEhgcLHS",     
		"AVEhgcHSG", 
	    "AVEhgcSCL", 
	    "AVEhgcCLGmore", 
	    
		"everAVEhgcLHS", 
		"everAVEhgcHSG", 
		"everAVEhgcSCL", 
		"everAVEhgcCLGmore", 
	  
		"AVEhgcHSGmore", 
		"everAVEhgcHSGmore", 

		    
		"everAVEdqageFifteen", 
		"everAVEdqageTwenty", 
		"everAVEdqageTwentyfive", 
		"everAVEdqageThirty", 
		

	    "everAVEadd",        
	    "everAVEqstock",   
	    "everAVEdq",   

	       "everAVEhgc",  		
	       "everAVEearnings", 
	       "everAVEc", 
	       "everAVEnetworth",  
              		   		
	       "everProbhgcHSG",                      
	       "everProbhgcSCL", 
	       "everProbhgcCLG", 
        
		   "everAVEeverdq",   
		   "everAVEdqave",   

	       "everAVEskill", 
	       "everAVEinitvalue", 
	       "everAVEvalueThirty", 
	       

   	    "AGEdqAVEadd",        
   	    "AGEdqAVEqstock",   
   	    "AGEdqAVEdq",   

   	       "AGEdqAVEhgc",  		
   	       "AGEdqAVEearnings", 
   	       "AGEdqAVEc", 
   	       "AGEdqAVEnetworth",  
              		   		
   	       "AGEdqProbhgcHSG",                      
   	       "AGEdqProbhgcSCL", 
   	       "AGEdqProbhgcCLG", 
        
   		   "AGEdqAVEAGEdqdq",   
   		   "AGEdqAVEdqave",   

   	       "AGEdqAVEskill", 
   	       "AGEdqAVEinitvalue", 
   	       "AGEdqAVEvalueThirty", 

	
	
	"AVEdqeverageTwentyfive",
	
	
	"AVEdqageFifteen", 
	"AVEdqageTwenty", 
	"AVEdqageTwentyfive", 
	"AVEdqageThirty", 

	"AVEhgcageFifteen", 
	"AVEhgcageTwenty", 
	"AVEhgcageTwentyfive", 
	"AVEhgcageThirty", 
	

	"AVEcageFifteen", 
	"AVEcageTwenty", 
	"AVEcageTwentyfive", 
	"AVEcageThirty", 


	"AVEnetworthageFifteen", 
	"AVEnetworthageTwenty", 
	"AVEnetworthageTwentyfive", 
	"AVEnetworthageThirty", 

	"AVEearningsageFifteen", 
	"AVEearningsageTwenty", 
	"AVEearningsageTwentyfive", 
	"AVEearningsageThirty", 


	"AVEvalueageFifteen", 
	"AVEvalueageTwenty", 
	"AVEvalueageTwentyfive", 
	"AVEvalueageThirty", 


	"AVEqstockageFifteen", 
	"AVEqstockageTwenty", 
	"AVEqstockageTwentyfive", 
	"AVEqstockageThirty", 
	
	
	"AVEaddageFifteen", 
	"AVEaddageTwenty", 
	"AVEaddageTwentyfive", 
	"AVEaddageThirty" 
		
                  
   	};
	
	int nstring = sizeof(vec_str)/sizeof(string); 
    const int nstat = nstring; 
	    
    double mat_stat[nstat][Num_cf]; 
    
    
    for(int i_n=0; i_n < Num_cf; i_n++){   
	
		
	int istat = 0; 

	  mat_stat[istat][i_n] = AVEagedeath_a15 [i_n];  istat++;           
	  mat_stat[istat][i_n] = AVEyearsleft_a15[i_n];  istat++;     
	  
	  mat_stat[istat][i_n] = everAVEagedeath_a15 [i_n];  istat++;           
	  mat_stat[istat][i_n] = everAVEyearsleft_a15[i_n];  istat++;  	        
      
      mat_stat[istat][i_n] = 100.0*(AVEhgc1315 [i_n] +  AVEhgc16more [i_n])     ; 			istat++;           
      mat_stat[istat][i_n] = 100.0*(AVEhgc1315_everdq[i_n] +  AVEhgc16more_everdq[i_n])    ; 			istat++; 

      mat_stat[istat][i_n] = 100.0*AVEhgc11less[i_n] ; 			istat++;     
      mat_stat[istat][i_n] = 100.0*AVEhgc12 [i_n]    ; 			istat++; 
      mat_stat[istat][i_n] = 100.0*AVEhgc1315 [i_n]  ; 			istat++;          
      mat_stat[istat][i_n] = 100.0*AVEhgc16more [i_n]; 			istat++; 
      
      mat_stat[istat][i_n] = 100.0*AVEhgc11less_everdq [i_n]; 			istat++;      
      mat_stat[istat][i_n] = 100.0*AVEhgc12_everdq [i_n]    ; 			istat++; 
      mat_stat[istat][i_n] = 100.0*AVEhgc1315_everdq  [i_n] ; 			istat++;          
      mat_stat[istat][i_n] = 100.0*AVEhgc16more_everdq[i_n] ; 			istat++;  

      mat_stat[istat][i_n] = 100.0*AVEhgc12more [i_n]       ; 			istat++;           
      mat_stat[istat][i_n] = 100.0*AVEhgc12more_everdq [i_n]; 			istat++; 

      
	{
	
	mat_stat[istat][i_n] = 100.0*everAVEdq_a15[i_n]; 			istat++; 
	mat_stat[istat][i_n] = 100.0*everAVEdq_a20[i_n]; 			istat++; 
	mat_stat[istat][i_n] = 100.0*everAVEdq_a25[i_n]; 			istat++; 
	mat_stat[istat][i_n] = 100.0*everAVEdq_a30[i_n]; 			istat++; 

	
	mat_stat[istat][i_n] = everAVEadd[i_n]; 			istat++; 
	mat_stat[istat][i_n] = everAVEqstock[i_n]; 			istat++; 
	mat_stat[istat][i_n] = 100.0*everAVEdq[i_n]; 			istat++; 		

	mat_stat[istat][i_n] = everAVEhgc[i_n]; 			istat++; 	
	mat_stat[istat][i_n] = everAVEearnings[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = everAVEc[i_n]/1.0; 				istat++; 
	mat_stat[istat][i_n] = everAVEnetworth[i_n]/1.0; 			istat++; 

	mat_stat[istat][i_n] = 100.0*everProbhgc12_a19[i_n]; 		istat++; 			
	mat_stat[istat][i_n] = 100.0*everProbhgc13_a21[i_n]; 		istat++; 
	mat_stat[istat][i_n] = 100.0*everProbhgc16_a25[i_n]; 		istat++; 

	mat_stat[istat][i_n] = 100.0*everAVEeverdq[i_n]; 			istat++; 	
	mat_stat[istat][i_n] = 100.0*everAVEdq1630[i_n]; 			istat++; 	
	
	mat_stat[istat][i_n] = everAVEskill[i_n]; 		    istat++; 			
	mat_stat[istat][i_n] = everAVEinitvalue[i_n]; 			istat++; 
	mat_stat[istat][i_n] = everAVEvaluea30[i_n]; 		    istat++; 		
		
	}


	{
	
	mat_stat[istat][i_n] = dqa30AVEadd[i_n]; 			istat++; 
	mat_stat[istat][i_n] = dqa30AVEqstock[i_n]; 			istat++; 
	mat_stat[istat][i_n] = 100.0*dqa30AVEdq[i_n]; 			istat++; 		

	mat_stat[istat][i_n] = dqa30AVEhgc[i_n]; 			istat++; 	
	mat_stat[istat][i_n] = dqa30AVEearnings[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = dqa30AVEc[i_n]/1.0; 				istat++; 
	mat_stat[istat][i_n] = dqa30AVEnetworth[i_n]/1.0; 			istat++; 

	mat_stat[istat][i_n] = 100.0*dqa30Probhgc12_a19[i_n]; 		istat++; 			
	mat_stat[istat][i_n] = 100.0*dqa30Probhgc13_a21[i_n]; 		istat++; 
	mat_stat[istat][i_n] = 100.0*dqa30Probhgc16_a25[i_n]; 		istat++; 

	mat_stat[istat][i_n] = 100.0*dqa30AVEdqa30dq[i_n]; 			istat++; 	
	mat_stat[istat][i_n] = 100.0*dqa30AVEdq1630[i_n]; 			istat++; 	
	
	mat_stat[istat][i_n] = dqa30AVEskill[i_n]; 		    istat++; 			
	mat_stat[istat][i_n] = dqa30AVEinitvalue[i_n]; 			istat++; 
	mat_stat[istat][i_n] = dqa30AVEvaluea30[i_n]; 		    istat++; 		
		
	}

	
	mat_stat[istat][i_n] = 100.0*AVEdqever_a25[i_n]; 		istat++; 
	
	
	mat_stat[istat][i_n] = 100.0*AVEdq_a15[i_n]; 		istat++; 
	mat_stat[istat][i_n] = 100.0*AVEdq_a20[i_n]; 		istat++; 
	mat_stat[istat][i_n] = 100.0*AVEdq_a25[i_n]; 		istat++; 
	mat_stat[istat][i_n] = 100.0*AVEdq_a30[i_n]; 		istat++; 

	mat_stat[istat][i_n] = AVEhgc_a15[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEhgc_a20[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEhgc_a25[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEhgc_a30[i_n]; 		istat++; 
	

	mat_stat[istat][i_n] = AVEc_a15[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEc_a20[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEc_a25[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEc_a30[i_n]/1.0; 		istat++; 


	mat_stat[istat][i_n] = AVEnetworth_a15[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEnetworth_a20[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEnetworth_a25[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEnetworth_a30[i_n]/1.0; 		istat++; 

	mat_stat[istat][i_n] = AVEearnings_a15[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEearnings_a20[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEearnings_a25[i_n]/1.0; 		istat++; 
	mat_stat[istat][i_n] = AVEearnings_a30[i_n]/1.0; 		istat++; 


	mat_stat[istat][i_n] = AVEvalue_a15[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEvalue_a20[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEvalue_a25[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEvalue_a30[i_n]; 		istat++; 


	mat_stat[istat][i_n] = AVEqstock_a15[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEqstock_a20[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEqstock_a25[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEqstock_a30[i_n]; 		istat++; 
	
	
	mat_stat[istat][i_n] = AVEadd_a15[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEadd_a20[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEadd_a25[i_n]; 		istat++; 
	mat_stat[istat][i_n] = AVEadd_a30[i_n]; 		istat++; 
		
   
       
    if (istat != nstring) { printf("\n Errors in write_all.cpp: write_input(): istat %3d != nstring %3d \n exiting ...\n", istat, nstring); exit(1); }
        
	}  // i_n 
		
	
	} // parafit != 1
	
	return 0;
}



int outfiles_all(string fdir, int i_n){
	
	get_stat(i_n); 
	
	string str; 
	char charname2[1024]; 
	
	str = fdir + "/moments_vec.txt"; 
	strcpy(charname2, str.c_str()); 
	printvec(charname2, SMM_moments, NUM_ismm);		

	str = fdir + "/estpar.txt"; 
	strcpy(charname2, str.c_str()); 
	printpar(charname2);  	
	
	
	str = fdir + "/educ.txt"; 
	strcpy(charname2, str.c_str()); 	
	FILE * feduc = fopen (charname2, "w"); 
	fprintf(feduc, "%.4f \n %.4f \n %.4f \n %.4f \n", AVEhgc12more[i_n], AVEhgc1315[i_n]+AVEhgc16more[i_n], AVEhgc16more[i_n], AVEhgc[i_n]); 
	fprintf(feduc, "%.4f \n %.4f \n %.4f \n %.4f \n", AVEhgc12more_everdq[i_n], AVEhgc1315_everdq[i_n]+AVEhgc16more_everdq[i_n], AVEhgc16more_everdq[i_n], everAVEhgc[i_n]); 
	fclose(feduc); 
	
	
	// The complete data (sim_model.txt) contains restricted access inforamtion and cannot be provided
	// Instead, we provide de-identified part of the data: sim_model_part1 for the purpose of replicating the main results of the paper 
	str = fdir + "/sim_model.txt"; 
	strcpy(charname2, str.c_str()); 
	write_smOpt(charname2, 0, 30-INIT_age); 	
	
	return 0; 

}



int CF_model (){
 	
 	
 	cc_taxPolicy = 0; 
    cc_taxPolicy_incomeTest = 1.0e+30; 	
	TAXrate_y  = 0.0;  
	TAXrate_dq = 0.0; 
	cbc_sub = 0.0;  
	gtr_sub   = 0.0;  
	gtr_csub  = 0.0; 
 	
 	   
  	ifstream infile; 	
	gsl_vector * vecx = gsl_vector_alloc(Nvar);		

	int i_n = 0; double dis; 

/*	
	{
		

	// *** fitted model ***	
	printf("\n\n *** Fitted Model *** \n\n");	
	i_n = 0;

    init_data();     
    init_sm0(0, NUM_smmPerson-1);
	infile.open("estpar_loglike_iter.txt"); readpar_loglike(500, infile);	infile.close(); 	
	infile.open("estpar_iter.txt"); 	    readpar(500, infile);	infile.close(); 
    get_coefEmaxstate(); 	// must have to make sure get_coefEmaxstate() is updated ... 
    
	
	time_t start0,end0;	double timediff; 
	time (&start0);
	
	getvec(vecx); getpar(vecx);
	dis = smm_obj(vecx,vecx); 
	
	time (&end0);	timediff = difftime (end0,start0);
	printf ("\n\n It took %f minutes to evaluate moment conditions, Initial distance = %f \n\n", timediff/60, dis);

	i_n = 0; 
    printf("\t\t i_n %d, dis %f ", i_n, dis);	
   
    outfiles_all("out_files", i_n); 

	i_n ++; 
	
    write_input(1);
	
	}
*/

	{
	
	int cf_ia = 0; 		
	
		i_n = 5; 
		cf_ia = 20 - INIT_age; 
		printf("\n\n *** %d: Effects of Education at age %d (make sure run benchmark model first) \n\n", i_n, cf_ia+INIT_age); 
	
		printf("\n check 1: mean hgc %d %f", cf_ia+INIT_age, my_mean(NUM_smmPerson, sm_educ [cf_ia])); 
				
		for (int iper = 0; iper <= NUM_smmPerson-1; iper++){
			int ia = cf_ia; 
			sm_educ [ia][iper] += 1; 	if ( sm_educ[ia][iper] < MIN_educ ) sm_educ[ia][iper] = MIN_educ; 
			if (sm_educ [ia][iper] > MAX_educ)	 sm_educ [ia][iper] = MAX_educ; 	
		}
		printf("\n check 1: mean hgc %d %f", cf_ia+INIT_age, my_mean(NUM_smmPerson, sm_educ [cf_ia])); 
		
	int MAX_initeduc_old = MAX_initeduc ;  
	MAX_initeduc = MAX_initeduc_old +1; 
	get_coefEmaxstate(); 
	SolveModel_EGM();	
						
	get_smOpt_EGM(0, NUM_smmPerson-1, cf_ia, NUM_age-1); 
	
		printf("\n check 1: mean hgc %d %f", cf_ia+INIT_age, my_mean(NUM_smmPerson, sm_educ [cf_ia])); 
			
    dis = get_sumf2();
    printf("\t\t i_n %d, dis %f ", i_n, dis);	
    outfiles_all("out_files_CF_EducPlus_3", i_n); 
	i_n ++; 
	
	
	MAX_initeduc = MAX_initeduc_old; 
	get_coefEmaxstate(); 
			
	}
	
	
	


	if (i_n > Num_cf){printf("\n\n Errors: i_n %d > Num_cf %d \n exiting ... \n", i_n, Num_cf); exit(1);}

	gsl_vector_free(vecx); 	

	
	write_input(0);  
	
	
	return 0; 
	
}


/*****************************************************************************************/

int prog_taxPolicy_model (int i_model, const int Niter_tax){
	
	if (i_model != 0 && i_model != 1 && i_model != 2) { printf("\n Error, unknown government transfer type \n exiting...\n"); exit(1); }
	
	// i_model = 0: subsidizing schooling 
	// i_model = 1: subsidizing at risk population  
	// i_model = 2: subsidizing non-smokers 
	
	double gtr_sub_new;
		
	// iteration for revenue neutral transfers
	int iter_tax = 0; 	int converge = 0; 
		
	while ( converge == 0 && iter_tax <= Niter_tax){
			
		iter_tax ++; 
		
		
		SolveModel_EGM();	
		
		
		// ... inner loop ... get_smOpt is fast while get_coefEmax is very time consuming 
		int Niter2_tax = 1 + 4 * (Niter_tax > 0);
		
		
		for (int iter2_tax = 0; iter2_tax < Niter2_tax; iter2_tax++) { 
		
			get_smOpt_EGM(0, NUM_smmPerson-1, 0, NUM_age-1);  // NUM_age-1
		
			
			gtr_sub_new = 0.0; 
			int n_sub = 0; // int n_de = 0; 
			double pop = 0.0; double pop_taxed = 0.0; 
			
			for (int ia = 0; ia <= NUM_age-1; ia++){
			for (int iper =0; iper < NUM_smmPerson; iper++){
				double cost_dq = sm_cprice[ia][iper]; 
				
				// population weight to adjust for population size due to death 
				double popwgt = 1.0; 
				if (ia >=1 ) popwgt = sm_probs[ia-1][iper]; 
				
				gtr_sub_new += popwgt * sm_dq[ia][iper] * cost_dq * TAXrate_dq; 
				
				double theta[] = {sm_theta[0][iper], sm_theta[1][iper]}; 
				
				int isub         = getGovSubRule(i_model, sm_de[ia][iper], sm_dq[ia][iper], theta, sm_educ[ia][iper], ia); 
				
				n_sub            += popwgt * isub; 
				
				pop       += popwgt; 
				pop_taxed += popwgt * sm_dq[ia][iper];
			}
			}		
			
			if (n_sub >= 1 ){
				gtr_sub_new /= (1.0*n_sub); 
			} else {
				gtr_sub_new = 0.0; // no one is eligible for the subsidy 
			}
			
			
			double maxv = 1.0; 
			double diff = 0.0; 
			
			diff  = fabs(gtr_sub_new - (cbc_sub+gtr_sub+gtr_csub) ); 
						
			maxv = fabs(cbc_sub+gtr_sub+gtr_csub);
			if (maxv < 1.0) maxv = 1.0; 			
			
			if (iter2_tax == 0)  printf("\n TAXrate_y %f, TAXrate_dq %f, \t iter_tax %2d, iter2_tax %d, relative diff %f: diff %f, gtr_sub_new %f, gtr_sub_old %f , pct subsidy %.2f  \t age 30 hgc %f, dq %.2f \n", TAXrate_y, TAXrate_dq, iter_tax, iter2_tax, diff/maxv, diff, gtr_sub_new, (cbc_sub+gtr_sub+gtr_csub), 1.0*n_sub/(1.0*NUM_age*NUM_smmPerson), my_mean(NUM_smmPerson, sm_educ[30 - INIT_age]), my_mean(NUM_smmPerson, sm_dq[30 - INIT_age])); // , my_mean(NUM_smmPerson, sm_cprice[ia])
							
		// check convergence:   
		if ( diff/maxv < 0.001 || diff < 10.0 ){
			
			if (iter_tax > 1 && iter2_tax == 0){
				
				converge = 1; 

                printf("\n\t TAXrate_dq %.2f, \t iter_tax %2d, iter2_tax %d, relative diff %f: diff %f, gtr_sub_new %f, gtr_sub_old %f , pct subsidy %.2f, pct taxed %.2f  \t age 30 hgc %f, dq %.2f \n", TAXrate_dq, iter_tax, iter2_tax, diff/maxv, diff, gtr_sub_new, (cbc_sub+gtr_sub+gtr_csub), 1.0*n_sub/pop, pop_taxed/pop, my_mean(NUM_smmPerson, sm_educ[30 - INIT_age]), my_mean(NUM_smmPerson, sm_dq[30 - INIT_age])); // , my_mean(NUM_smmPerson, sm_cprice[ia])
				
			}
			
			break; // break the loop of iter2_tax 
			
		} else { 
			
			
			// ... update ... 
				if (i_model == 0){
					cbc_sub = gtr_sub_new; 
					gtr_sub = 0.0; 
					gtr_csub = 0.0; 
				} else if (i_model == 1){
					cbc_sub  = 0.0;  				
					gtr_sub  = gtr_sub_new; 
					gtr_csub = 0.0; 
				} else {
					cbc_sub  = 0.0; 		
					gtr_sub  = 0.0;  
					gtr_csub = gtr_sub_new;
				} 
				
		} // check for convergence 
				
		
		} // iter2_tax
							
	} 
	
	if (converge != 1 && Niter_tax > 0){ printf("\n\n Warnings: gov subsidy does not converge !!!!\n\n"); }
	
	return 0; 
	
}



int get_optTax_model (const string fdir, const int i_n, const int i_model, double maxtaxq, int ngrid){

		// ... read education from noDq simulation
		ifstream infile; 
	    infile.open("out_files_CFnoDq/educ.txt"); 	
	    if(!infile) {cerr<<"\n Error reading the file out_files_CFnoDq/educ.txt \n"; exit(1); }
		double target_educ; double educ_CFnoDq[8]; 
		for (int iv = 0; iv < 8; iv++){
		infile >> educ_CFnoDq[iv];  // college attendance age 21, 4-year college grad age 25, state: hgc age 30
		}
		infile.close(); 
		
		int dummy_everdq[NUM_smmPerson]; 
		// ... read dummy_everdq from benchmark model 
	    infile.open("out_files/dummy_everdq.txt"); 	
	    if(!infile) {cerr<<"\n Error reading the file out_files/dummy_everdq.txt \n"; exit(1); }
		for (int iper = 0; iper < NUM_smmPerson; iper++){
		infile >> dummy_everdq[iper]; 
		}		
		infile.close(); 
	
		target_educ = educ_CFnoDq[5]; // hgc >= 13

		
	printf("\n target_educ %f , ", target_educ);

		
	
	
	gsl_vector * vecx = gsl_vector_alloc(Nvar);	


	
	cc_taxPolicy = 1; 
    cc_taxPolicy_incomeTest = 1.0e+30; 	
	
	const int Niter_tax = 20;     
	
	
	{
		
	string str; 
	char charname2[1024]; 
	
	FILE * outfile1; FILE * outfile2; 	
	
	str = fdir + "/ALL_outcomes.txt"; 	strcpy(charname2, str.c_str()); 
	outfile1         = fopen (charname2, "w");
	str = fdir + "/OPT_outcomes.txt"; 	strcpy(charname2, str.c_str()); 	
	outfile2         = fopen (charname2, "w");
	

    init_data();     
    init_sm0(0, NUM_smmPerson-1);
	infile.open("estpar_loglike_iter.txt"); readpar_loglike(500, infile);	infile.close(); 	
	infile.open("estpar_iter.txt"); 	    readpar(500, infile);	infile.close(); 
    get_coefEmaxstate(); 	// must have to make sure get_coefEmaxstate() is updated ... 
	
	getvec(vecx); getpar(vecx);

	
	double opt_tax_dq = 0.0;
	double opt_tax_y = 0.0; 	
	double opt_cbc_sub, opt_gtr_sub, opt_gtr_csub; 
	double opt_dqa30 = 0.0; double opt_hgc30 = 0.0; double opt_hgc12more = 0.0; double opt_hgc16more = 0.0; 
	double opt_dqa30_everdq, opt_hgc30_everdq, opt_hgc12more_everdq, opt_hgc13more_everdq, opt_hgc16more_everdq;

	double diff_hgc30 = 1.0e+30; 

	
	double dis;

	
    printf("\n\n *** %d: Optimal Corrective Taxation, tuition subsidy i_model %d \n\n", i_n, i_model);

	cbc_sub = 0.0;  
	gtr_sub = 0.0; 
	gtr_csub = 0.0; 
	
	TAXrate_y  = 0.0;  
	
	double diff_hgc30_new = 10000.0; 

	
	for (int igrid=1; igrid<= ngrid; igrid++)
    {
		
			
		TAXrate_dq = igrid * maxtaxq / ngrid; 
        
        printf("\n\n igrid %d, taxrate_dq %f", igrid, TAXrate_dq);
        
		// find equilibrium cash transfer for each age
			prog_taxPolicy_model(i_model, Niter_tax); 
			
            // sm_educcat[i2530][4][iper] is calculated in get_sumf2 
            int i2530 = 1; 
            int iage30 = 30 - INIT_age; 
			
			const double hgc30 = my_mean(NUM_smmPerson, sm_educ[iage30]); 	
			const double dqa30        = my_mean(NUM_smmPerson, sm_dq[iage30]);  		
            
			dis = get_sumf2();     
						
            double hgc11less = my_mean(NUM_smmPerson, sm_educcat[i2530][0]);            
            double hgc12     = my_mean(NUM_smmPerson, sm_educcat[i2530][1]); 
            double hgc1315   = my_mean(NUM_smmPerson, sm_educcat[i2530][2]);             
            double hgc16more = my_mean(NUM_smmPerson, sm_educcat[i2530][3]);

			double dqa30_everdq = my_mean_conditional(NUM_smmPerson, sm_dq[iage30],   dummy_everdq, 1,  1);  			
			double hgc30_everdq = my_mean_conditional(NUM_smmPerson, sm_educ[iage30], dummy_everdq, 1,  1); 
			
            double hgc11less_everdq = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][0], dummy_everdq, 1,  1);        
            double hgc12_everdq     = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][1], dummy_everdq, 1,  1); 
            double hgc1315_everdq   = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][2], dummy_everdq, 1,  1);           
            double hgc16more_everdq = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][3], dummy_everdq, 1,  1); 
		
		diff_hgc30_new = ( target_educ - (hgc1315_everdq+hgc16more_everdq)); 	
		
		
		if ( fabs (diff_hgc30_new ) < diff_hgc30  )
		{
            
            diff_hgc30 = diff_hgc30_new;
            
			opt_tax_y = TAXrate_y; 
			opt_tax_dq = TAXrate_dq;
			opt_cbc_sub = cbc_sub; 
			opt_gtr_sub = gtr_sub; 
			opt_gtr_csub = gtr_csub; 
			
			opt_dqa30 = dqa30; 
			opt_hgc30 = hgc30;
			opt_hgc12more = 1.0 - hgc11less; 
			opt_hgc16more = hgc16more; 

			opt_dqa30_everdq = dqa30_everdq; 
			opt_hgc30_everdq = hgc30_everdq;
			opt_hgc12more_everdq = 1.0 - hgc11less_everdq; 
			opt_hgc16more_everdq = hgc16more_everdq;
            opt_hgc13more_everdq =hgc1315_everdq +hgc16more_everdq;
         }
		

		 	
			printf("\t i_n %d, dis %f, tax_dq %f \t diff_hgc %f, %f \t Ever smokers: HSG+ %.4f, SC+ %.4f, CG+ %.4f \n\n", i_n, dis, TAXrate_dq, diff_hgc30_new, diff_hgc30, 1.0-hgc11less_everdq, hgc1315_everdq+hgc16more_everdq, hgc16more_everdq  );
			
			// overall population
			fprintf(outfile1, "%d \t %d \t %.8f \t %.8f \t %.8f \t %.8f \t ", 0, igrid, TAXrate_dq, cbc_sub, gtr_sub, gtr_csub);	
			fprintf(outfile1, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \t %.8f  \t %.8f \n",  dqa30, hgc30, 1.0-hgc11less, hgc11less, hgc12, hgc1315, hgc16more );				
			// among ever smokers
			fprintf(outfile1, "%d \t %d \t %.8f \t %.8f \t %.8f \t %.8f \t ", 1, igrid, TAXrate_dq, cbc_sub, gtr_sub, gtr_csub);	
			fprintf(outfile1, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \t %.8f  \t %.8f \n",  dqa30_everdq, hgc30_everdq, 1.0-hgc11less_everdq, hgc11less_everdq, hgc12_everdq, hgc1315_everdq, hgc16more_everdq );				
		
		
	}
	fclose(outfile1); 
	
	
		
	TAXrate_y  = opt_tax_y ;
	TAXrate_dq = opt_tax_dq; 
	cbc_sub    = opt_cbc_sub ; 
	gtr_sub    = opt_gtr_sub ; 
	gtr_csub   = opt_gtr_csub ; 
			
	dis = smm_obj(vecx,vecx); 	
    
    printf("\t\t i_n %d, dis %f, hgc %.2f, hgc13more %f, target %.2f, taxrate %f ", i_n, dis, my_mean(NUM_smmPerson, sm_educ[30 - INIT_age]), opt_hgc13more_everdq, target_educ, TAXrate_dq);

   			 outfiles_all(fdir, i_n);
		
	        // sm_educcat[i2530][4][iper] is calculated in get_sumf2 
            int i2530 = 1; 
            int iage30 = 30 - INIT_age; 
		
			const double hgc30 = my_mean(NUM_smmPerson, sm_educ[iage30]); 	
			const double dqa30        = my_mean(NUM_smmPerson, sm_dq[iage30]);    
						
            double hgc11less = my_mean(NUM_smmPerson, sm_educcat[i2530][0]);            
            double hgc12     = my_mean(NUM_smmPerson, sm_educcat[i2530][1]); 
            double hgc1315   = my_mean(NUM_smmPerson, sm_educcat[i2530][2]);             
            double hgc16more = my_mean(NUM_smmPerson, sm_educcat[i2530][3]);

			double dqa30_everdq = my_mean_conditional(NUM_smmPerson, sm_dq[iage30],   dummy_everdq, 1,  1);  			
			double hgc30_everdq = my_mean_conditional(NUM_smmPerson, sm_educ[iage30], dummy_everdq, 1,  1); 
			
            double hgc11less_everdq = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][0], dummy_everdq, 1,  1);        
            double hgc12_everdq     = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][1], dummy_everdq, 1,  1); 
            double hgc1315_everdq   = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][2], dummy_everdq, 1,  1);           
            double hgc16more_everdq = my_mean_conditional(NUM_smmPerson, sm_educcat[i2530][3], dummy_everdq, 1,  1); 
			            
			fprintf(outfile2, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \n ", TAXrate_y, TAXrate_dq, cbc_sub, gtr_sub, gtr_csub);	
			fprintf(outfile2, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \t %.8f  \t %.8f \n",  dqa30, hgc30, 1.0-hgc11less, hgc11less, hgc12, hgc1315, hgc16more );				
			fprintf(outfile2, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \t %.8f  \t %.8f \n",  dqa30_everdq, hgc30_everdq, 1.0-hgc11less_everdq, hgc11less_everdq, hgc12_everdq, hgc1315_everdq, hgc16more_everdq );				
			
			int igrid = floor(TAXrate_dq/0.05); 
			// overall population
			fprintf(outfile2, "%d \t %d \t %.8f \t %.8f \t %.8f \t %.8f \t ", 0, igrid, TAXrate_dq, cbc_sub, gtr_sub, gtr_csub);	
			fprintf(outfile2, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \t %.8f  \t %.8f \n",  dqa30, hgc30, 1.0-hgc11less, hgc11less, hgc12, hgc1315, hgc16more );				
			// among ever smokers
			fprintf(outfile2, "%d \t %d \t %.8f \t %.8f \t %.8f \t %.8f \t ", 1, igrid, TAXrate_dq, cbc_sub, gtr_sub, gtr_csub);	
			fprintf(outfile2, "%.8f \t %.8f \t %.8f \t %.8f \t %.8f \t %.8f  \t %.8f \n",  dqa30_everdq, hgc30_everdq, 1.0-hgc11less_everdq, hgc11less_everdq, hgc12_everdq, hgc1315_everdq, hgc16more_everdq );				

			fclose(outfile2); 
	
	
	if (i_model == 0){ // cash transfer for schooling 
			
	FILE * input_tax; 
	input_tax  = fopen ("input_tax.tex", "w"); 
	double elasticity[3]; 
	{
        printf("\n dq [%f, %f, %f], benchmark [%f, %f, %f], tax %f \n", AVEdq_a20[i_n], AVEdq_a25[i_n], AVEdq_a30[i_n], AVEdq_a20[0], AVEdq_a25[0], AVEdq_a30[0], opt_tax_dq);
	elasticity[0] = ((AVEdq_a20[i_n] - AVEdq_a20[0])/AVEdq_a20[0])/opt_tax_dq;
	elasticity[1] = ((AVEdq_a25[i_n] - AVEdq_a25[0])/AVEdq_a25[0])/opt_tax_dq;	
	elasticity[2] = ((AVEdq_a30[i_n] - AVEdq_a30[0])/AVEdq_a30[0])/opt_tax_dq; 
	}
   	string str; 
	string vec_strNumb[] = { "CFZero", "CFOne", "CFTwo", "CFThree", "CFFour", "CFFive", "CFSix", "CFSeven", "CFEight"}; 	
	string vec_str[] = {"optTAXratey", "optTAXratedq", "optcbcsub", "optgtrsub", "optgtrCsub", "optdq", "opthgc", "opthgcHSGmore", "opthgcCLGmore", "optdqever", "opthgcever", "opthgcHSGmoreever", "opthgcCLGmoreever", "optElasticityTwenty", "optElasticityTwentyfive", "optElasticityThirty"    }; 
	double vec_tax[16]; 

	vec_tax[0] = opt_tax_y;  vec_tax[1] = opt_tax_dq; vec_tax[2] = opt_cbc_sub; vec_tax[3] = opt_gtr_sub; vec_tax[4] = opt_gtr_csub; 	
	vec_tax[5] = opt_dqa30;  vec_tax[6] = opt_hgc30;  vec_tax[7] = opt_hgc12more; vec_tax[8] = opt_hgc16more;
	vec_tax[9] = opt_dqa30_everdq;  vec_tax[10] = opt_hgc30_everdq;  vec_tax[11] = opt_hgc12more_everdq; vec_tax[12] = opt_hgc16more_everdq;
	vec_tax[13] = elasticity[0]; vec_tax[14] = elasticity[1]; vec_tax[15] = elasticity[2];
	
	for (int ivt=0; ivt < 16; ivt++){
	   	str = vec_str[ivt] + vec_strNumb[i_model]; 
//		fprintf(input_tax,  "	\\newcommand{	\\%s			}{  %.0f } \n",	str.c_str(), vec_tax[ivt] ); 
		fprintf(input_tax,  "	\\newcommand{	\\PCT%s			}{  %.0f } \n",	str.c_str(), vec_tax[ivt]*100.0 ); 
		fprintf(input_tax,  "	\\newcommand{	\\d%s			}{  %.1f } \n",	str.c_str(), vec_tax[ivt] ); 	
		fprintf(input_tax,  "	\\newcommand{	\\dd%s			}{  %.2f } \n",	str.c_str(), vec_tax[ivt] ); 							
		fprintf(input_tax,  "	\\newcommand{	\\ddd%s			}{  %.3f } \n",	str.c_str(), vec_tax[ivt] ); 							
		fprintf(input_tax,  "	\\newcommand{	\\dddd%s			}{  %.4f } \n",	str.c_str(), vec_tax[ivt] ); 							
	}
    fclose(input_tax); 

		}
	
	} 

	gsl_vector_free(vecx); 	
	
	
	return 0; 

}



int CF_optTax_model (){
	
	// i_model controls the type of subsidy 
	int i_n, i_model; 
	

	ifstream infile; 	
	gsl_vector * vecx = gsl_vector_alloc(Nvar);		
	double dis; 
	
	// *** No smoking ***
	i_n = 1; 
	printf("\n\n *** %d: No smoking behavior \n\n", i_n); 
	init_data(); 
	init_sm0(0, NUM_smmPerson-1); 
	infile.open("estpar_loglike_iter.txt"); readpar_loglike(500, infile);	infile.close(); 	
	infile.open("estpar_iter.txt"); 	    readpar(500, infile);	infile.close(); 
    get_coefEmaxstate(); 	// must have to make sure get_coefEmaxstate() is updated ... 

		
	cc_taxPolicy = 1; 
    cc_taxPolicy_incomeTest = 1.0e+30; 	
	TAXrate_y  = 0.0;  
	TAXrate_dq = 1.0e+30; 
	cbc_sub = 0.0;  
	gtr_sub   = 0.0;  
	gtr_csub  = 0.0; 


	
	getvec(vecx); 
	dis = smm_obj(vecx,vecx); 
	
	printf("\t\t i_n %d, dis %f ", i_n, dis);	
	outfiles_all("out_files_CFnoDq", i_n); 	
	printf("\n\t\t HSG+ %.4f, SC+ %.4f, CG+ %.4f, mean hgc %.4f \n", AVEhgc12more[i_n], AVEhgc1315[i_n]+AVEhgc16more[i_n], AVEhgc16more[i_n], AVEhgc[i_n]); 
	

	// initialization 	
	TAXrate_dq = 0.35; 
		
		
	
	write_input(0);  
	
	
	i_n = 2; 
    i_model = 0; // cash transfer for attending school 
	get_optTax_model ("out_files_optTax", i_n, i_model, 0.35, 7); 		
	

	
	i_n = 4; 
    i_model = 2; // cash transfer for non-smokers 
	get_optTax_model ("out_files_optTax2", i_n, i_model, TAXrate_dq, 1); 	
	
	
	i_n = 3; 
    i_model = 1; // cash transfer for at-risk population 
	get_optTax_model ("out_files_optTax1", i_n, i_model, TAXrate_dq, 1); 	



 	cc_taxPolicy = 0; 
    cc_taxPolicy_incomeTest = 1.0e+30; 	
	TAXrate_y  = 0.0;  
	TAXrate_dq = 0.0; 
	cbc_sub = 0.0;  
	gtr_sub   = 0.0;  
	gtr_csub  = 0.0; 



	write_input(0);  
	

	if (i_n >= Num_cf){printf("\n\n CF_optTax_model () Errors: i_n %d > Num_cf %d \n exiting ... \n", i_n, Num_cf); exit(1);}


	// *** fitted model ***	
	printf("\n\n *** Check Fitted Model *** \n\n");	
//	i_n = 0;

    init_data();     
    init_sm0(0, NUM_smmPerson-1);
	infile.open("estpar_loglike_iter.txt"); readpar_loglike(500, infile);	infile.close(); 	
	infile.open("estpar_iter.txt"); 	    readpar(500, infile);	infile.close(); 
    get_coefEmaxstate(); 	// must have to make sure get_coefEmaxstate() is updated ... 
    
	
	time_t start0,end0;	double timediff; 
	time (&start0);
	
	getvec(vecx); getpar(vecx);
	dis = smm_obj(vecx,vecx); 
	
	time (&end0);	timediff = difftime (end0,start0);
	printf ("\n\n It took %f minutes to evaluate moment conditions, Initial distance = %f \n\n", timediff/60, dis);


	gsl_vector_free(vecx); 
			
	return 0; 
	
}
